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Abstract 

Structural health monitoring (SHM) for the detection of damage in aerospace materials is an important 
area of research at NASA. Ultrasonic guided Lamb waves are a promising SHM damage detection technique 
since the waves can propagate long distances. For complicated flaw geometries experimental signals can be 
difficult to interpret. High performance computing can now handle full 3-dimensional (3D) simulations of 
elastic wave propagation in materials. We have developed and implemented parallel 3D elastodynamic finite 
integration technique (3D EFIT) code to investigate ultrasound scattering from flaws in materials. EFIT 
results have been compared to experimental data and the simulations provide unique insight into details of 
the wave behavior. This type of insight is useful for developing optimized experimental SHM techniques. 3D 
EFIT can also be expanded to model wave propagation and scattering in anisotropic composite materials. 
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1 Introduction 

Developing optimized flaw detection techniques for aerospace materials is an important area 
of research in the fields of nondestructive evaluation (NDE) and structural health monitor- 
ing (SHM). As more advanced materials, such as composites, are developed and utilized in 
the aerospace industry, efficient and accurate damage detection techniques will be necessary 
to ensure vehicle/component safety and functionality. A benefit of SHM on-board sensor 
techniques over traditional ground based NDE is the potential detection of in-flight material 
damaging events (which could be particularly critical for future long-duration spaceflight 
vehicles). Current techniques used for finding flaws in aerospace materials include visual 
inspection, shearography, acoustic emission, eddy current testing, ultrasonics, and other 
methods. A review of SHM techniques currently used to find damage in composite air- 
craft material can be found in a recent paper by Diamanti [1]. An advancement in current 
techniques will be necessary for the realization of economical, efficient, and accurate SHM 
systems for aerospace vehicles. 

Guided Lamb waves are particularly ideal for detecting flaws in the plate and pipe-like 
structures that make up aerospace vehicles because the waves can propagate over long dis- 
tances (tens of meters) and can be sensitive to some flaw types that are difficult to find using 
traditional ultrasound techniques [2, 3]. Additionally, the use of multiple Lamb wave modes 
could allow for flaw detection throughout the depth of a material since each mode has a dif- 
ferent depth dependent displacement curve [4], Active sensors such as embedded/attached 
piezoelectric wafers, interdigital transducers, and macrofiber composites may be ideal for 
generating and detecting Lamb waves in SHM applications, where light-weight (and poten- 
tially flexible) sensors are required [5-7]. Several review articles discuss various aspects of 
the use of Lamb waves in SHM [5, 8]. 



Lamb waves are created due to the stress-free boundaries of the plate-like material. A 
coupling of longitudinal and shear waves leads to the creation of any number of dispersive 
Lamb wave modes, which separate into antisymmetric and symmetric modes in an isotropic 
medium [4], Although Lamb waves are described using the same elastodynamic equations as 
bulk waves, the presence of the plate boundaries leads to dispersion relations that must be 
solved numerically. The use of Lamb waves for nondestructive testing can be found in the 
literature as early as the 1950s [9], yet the complicated behavior of these dispersive waves 
has resulted in an extensive ongoing body of literature on guided wave based flaw detection. 
Techniques that excite multiple Lamb wave modes are much simpler than those required to 
create only one or two desired modes [8]. However, the generation of multiple Lamb wave 
modes leads to complicated wave behavior that can require sophisticated signal processing 
techniques to make sense of the data [10]. 

One way to provide insight into the complex behavior of Lamb wave interaction with flaws 
is to use numerical methods to simulate the Lamb wave propagation and scattering. The 
recent progress in computational power and availability allows us to run full 3-dimensional 
(3D) Lamb wave simulations using a computing cluster. In this paper we discuss the devel- 
opment, verification, and validation of an isotropic 3D simulation code that can be used to 
investigate Lamb wave behavior and interaction with material damage. One example of a 
longer term goal for these types of simulations would be to provide insight into optimal ultra- 
sonic sensor configurations for a specific SHM application. We also discuss the extension of 
the simulation code to include anisotropic behavior for modeling Lamb waves in composites. 
In the following section we describe the specific numerical method, called finite integration 
technique (FIT), that we have used to implement 3D Lamb wave simulations. 

2 Finite Integration Technique 

FIT was developed by Weiland in 1977 for electrodynamics modeling, and was applied to 
elastodynamics by Fellinger in the early 1990s [11, 12]. Since that time, the elastodynamic 
finite integration technique (EFIT) has been used to explore elastic wave behavior in a 
variety of applications. Marklein’s review paper gives a thorough overview of the history 
of FIT [13]. EFIT has been used over the last couple decades to investigate nondestructive 
evaluation applications such as scattering from flaws in plates, pipes, and bulk material 
[14, 15]. The technique has also been implemented to explore elastic wave propagation in 
anisotropic heterogeneous media for applications ranging from wave scattering from surface 
breaking cracks in steel welds to wave propagation in composites [16-18]. However, the 
majority of prior publications relating to EFIT have been limited to 2D simulations. 

We chose FIT over other numerical techniques, such as finite difference and finite element 
method, for several reasons: 1) FIT naturally requires staggered spatial and temporal grids 
which leads to stability, 2) boundary conditions are easily incorporated into FIT, and 3) the 
mathematical analysis is straight-forward and leads to equations that are easy to implement 
in any programming language [13]. Full 3D simulations of elastic waves in real-world sized 
simulation spaces are extremely computationally demanding and require access to computa- 
tional resources such as large computing clusters. We used Message Passing Interface (MPI) 
to create a customized parallel EFIT code that runs efficiently on a computing cluster at 



NASA Langley. We used a ID virtual topology for parallelization, breaking up the simula- 
tion space in the x± direction. This partitioning means that all equations requiring stress 
or velocity values in the +/ — x\ direction must have those values passed across processor 
boundaries, see figure 1. 


Variable: 



Required by: Send direction: 

yRight 


* Ti 3 ,T 12 /Vj 


T n ,v 2 ,v 3 



Figure 1: Diagram showing the variable passing required for EFIT with a ID virtual topology. 


Although EFIT equations can be found in the literature, below we give a brief mathe- 
matical description of the homogeneous EFIT equations for the reader’s convenience. Inho- 
mogeneous EFIT equations (which were used for the model verification in section 2.1) are 
available in several papers, such as [19], and will not be repeated here due to their additional 
length. Elastic wave propagation in isotropic media is described by Cauchy’s equation for 
momentum and Hooke’s Law. Taking the time derivative of these equations (in Cartesian 
coordinates) leads to nine equations for the velocities, v, and normal and shear stresses, T*j. 
In FIT the equations are integrated over a control volume, and the integrals are approxi- 
mated in a discrete form (details of this process can be found in [13]). To correctly transfer 
the equations into a discretized space, we place the stress and velocity components onto a 
discrete grid in the manner described by Fellinger [12]. For inhomogeneous material the 
Lame parameters p and A, and the density, p, are also discretized leading to equations that 
can be found in the papers mentioned above. Below we list the nine discretized homoge- 
neous EFIT equations describing wave propagation in isotropic material, where we assume 
Aaq = Ax 2 = Ax 3 [13]: 


■An) 


An) 


A n ) 


(t) 

(t) 

(t) 


pAx - 
1 


pAx - 
1 

pAx 


T^ +il \t) - T^\t) + T$(t) - T^~ i 2 \t) + T$(t) - T^~ £ 3 \t) 
+ T% +& 2 \t) 


FA) 


-,(71 — Xl) 
12 


T)s’(t) - Ttr"(t) + T"\t) - T . *”-“'(«) + r' 3 ”™'(«) - Tg'(t) + f 3 , 


(n+x 2) 
22 

(n-x i)/,\ , rp(n) 




+ /i 


+/ 2 , 


( 1 ) 


(n-x 2 ) , 


-,(n+x 3 ) , 


A n ), 


ifr’W 

T%\t) 

T£'(t) 


Ax 


(A^ + 2^)(v^\t) - vf-^it)) + X {n \vi n \t) - vt~* 2 \t) + vi n \t) 


(n-x 3 ) 


(*)) 


Ax 


(A (n) + 2^)(vi n \t) - v { 2 n ~ i2 \t)) + X^(v^(t) - v { r £l \t) + vi n \t) 


—u 


(n-x 3) 


(*)) 


( 2 ) 


Ax L 


(A^ + 2^)(v { 3 n \t) - v ( 3 n ~ i3 \t)) + X^(v^\t) - v^- £l \t) + vi n \t) 


(n-x 2) 


(*)) 



_J1_ 

/ (n+X 2 ) 

- v {n) ) 

+ 


(■ n ) 
— IV 

Ax 

\ U 1 


\ u 2 

v 2 


_J/_ 

Ax 

\v[ n+i3) 

(n)\ 

- v\ ’) 

+ 

(vi n+£l) 

( n ) 

-^3 

T'S'(t) 

_ d 

( (n+x 3 ) 

- V {n) ) 

+ 

( (n+x 2 ) 

- v {n) 

Ax 

\ u 2 

u 2 ) 

\ u 3 

u 3 


(3) 


where the notation (n + xi) represents one spatial step in the +fi direction beyond the 
current position, n, and the Lame parameters relate to the longitudinal, cl, and transverse, 
Ct- speeds of sound by cl = \/(A + 2 [i)/p and ct = The equations are discretized 

in time using central differences, which results in the velocity and stress components being 
staggered in time by (At/ 2) [20]: 
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where At and At/2 are full and half time steps. For each time step equations (1) - (3) are 
solved for all points in the simulation space. Using equations (1) through (5) the simulation 
proceeds forward in time. 

Stability conditions require minimum time and spatial step sizes. We will use the same 
size spatial step in all directions since the discrete space is broken up into cubes. The 
reported ideal number of points per wavelength varies in the literature, ranging from A/8 to 
A/15 [21, 22], The Courant-Friedrichs-Levy condition is used to determine the numerically 
stable time step size (with Aaq = Ax 2 = Ax 3 = Ax) [20, 22]: 
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where c ma x and c„ nn are the maximum and minimum wave speed, f max is the maximum 
frequency, and b is a constant with the criteria that b > 8. For the case of an elastic solid 
in air, stress-free boundary conditions must be applied at the surfaces. Equations related to 
stress-free boundaries can be found in [20]. 



2.1 Simulation Code Verification 

Implementation of the parallel 3D EFIT code was verified through various comparisons to 
analytical scattering solutions [23]. For example, we compared results from EFIT simulations 
of backscattering from an elastic sphere embedded in an elastic solid to a frequency domain 
analytical backscattering solution [24], The EFIT time domain data was transformed into 
frequency domain for the comparison. Figure 2 shows the comparison for a broadband 
incident wave scattering from an aluminum sphere embedded in brass, where the spatial 
step size was A/20. Note that since we did not include a perfectly matched layer at the 
simulation edges, the EFIT result differs slightly from the analytical result due to gating 
issues (i.e. edge reflections had to be gated out in time, thus also cutting out some of the 
backscattered waveform). We also point out that in order to minimize gating issues while 
maintaining a small step size, this comparison required a simulation space of 800 X 800 X 
800 steps (512 million matrix elements for each time step, totaling around 0.8 TB of data) 
and took approximately 50 hours of runtime on 32 processors (1600 CPU hours). 



Figure 2: Comparison of frequency domain EFIT (dotted line) and analytical solution (solid line) backscat- 
tering results for an aluminum sphere embedded in brass, where ka is wavenumber multiplied by sphere 
radius. 


3 Comparison to Experiment 

We compared EFIT simulation and experimental results for Lamb wave interaction with a 
void-type flaw in an aircraft grade aluminum plate. These simulation results were compared 
to experimental data for code validation purposes and to explore the type of insight that 
the 3D simulations can provide in relation to experimental results. The experimental setup 
consisted of a rounded-rectangle flat-bottom hole milled into a 305 mm x 305 mm x 3.154 
mm aluminum 2024 alloy plate with a rounded rectangle flaw of 76 mm in length and 30 
mm in width. Experimental data was taken with a circular 6.35 mm (0.25 inch) diameter 
transducer with a center frequency of 2.15 MHz. The generated incident wave was a 5-cycle 




sine wave tone burst. Glycerin was used as a couplant between a delay line and the plate 
surface. The transmitting and receiving transducers stepped through 100 locations spaced 
by 2 mm increments along both x and y axes of the plate, to cover a total scan area of 198 
mm x 198 mm. This pitch-catch setup was used to collect the experimental data for all 
corresponding ray paths as the flaw was incrementally milled deeper. Data was taken at 
sixteen void depths (including 0 and 100% material loss). 

For the corresponding EFIT simulations we incorporated the exact flaw shape and size. 
An incident wave of frequency 2.15 MHz and 5 cycles in length was introduced into the 
simulated plate. Each simulation took around 1400 CPU hours, and was completed in only 
44 wall-clock hours with the parallelized code running on 32 processors. An example of 
EFIT simulation results is shown in figure 3. The images show the out-of-plane motion 
only (plotting in-plane motion would also be useful for visualizing certain modes). Another 
technique we used for the comparisons was to track Lamb mode conversion in the simulation 
output using a short time Fourier transform (STFT) in spatial domain. The technique creates 
plots of wavenumber versus position, and more details can be found in [23]. In this paper 
we will only give an overview of our results from numerous comparisons between experiment 
and simulation output. 
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Figure 3: Example of EFIT simulation results (out-of-plane motion only) for a 3.154 mm thick aluminum 
plate with a 2.14 mm deep rounded rectangle flat bottom hole. The images show simulation output at a single 
point in time: a) 3D view showing the plate top and sides (flaw region outlined in red), b) 2D horizontal 
slice taken through the 3D space directly beneath the void showing waves propagating in the thinned region, 
c) 2D vertical slice showing the waves through the plate thickness (flaw region outlined in red). 


The simulations showed that for some transducer positions and void depths certain Lamb 
wave modes were significantly disrupted due to interference effects created by waves interact- 


ing with the flaw. In those cases the simulations helped explain the lack of detected modes 
in the experimental data. For other transducer locations we saw surprising results in the 
simulations which helped explain features in experimental data. Figure 4 shows simulation 
results, a spatial STFT plot, and experimental data for a case where the EFIT simulations 
predicted the creation of a high amplitude AO/SO Lamb mode after the waves propagated 
beneath the void. Simulation images showed that this high amplitude mode was most likely 
due to constructive interference created by the exact shape and size of the flaw region. As 
shown in the figure, a high amplitude mode does indeed appear in the experimental data at 
the time expected from simulation results. Note that all locations of the Lamb mode labels 
in the figure showing experimental data were chosen based only on EFIT simulation results. 
The EFIT predicted mode arrival times for this example, and for many other cases we in- 
vestigated, appear to match well with experimental data. Results for several comparisons 
between EFIT and experimental results can be found in [23]. The isotropic EFIT code will 
be used in the future to investigate other flaw types of interest. In the following section 
we briefly describe the extension of EFIT code to include anisotropic behavior for modeling 
wave behavior in composites. 
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Figure 4: a) EFIT result showing a 2D slice beneath a 2.68 mm deep flat-bottom hole at a single point in 
time, b) Wavenumber versus position plot created from EFIT simulation data, corresponding to the A-scan 
position shown by red line in a) and to that same point in time. The region beneath the void falls within 
the white rectangle and Lamb wave modes emerging post-flaw are labeled, c) Denoised experimental data 
for the unflawed plate and the plate with a 2.68 mm deep flat-bottom hole. The Lamb wave mode label 
locations in c) were chosen based on the wavenumber versus position plots created using EFIT simulation 
data. 




4 Anisotropic EFIT 


As mentioned earlier, anisotropic versions of EFIT have been implemented previously by 
a few authors [16-18]. Most of the reported simulation results, however, were limited to 
2D. Access to a large computing cluster will allow us to run 3D anisotropic elastic waves 
simulations. Incorporating anisotropy leads to significantly more complicated discretized 
stress equations compared to the homogeneous equations listed in section 2. The added 
complexity results from the fact that for a general triclinic anisotropic material the elastic 
matrix contains 21 different elastic constants [4], To demonstrate this increase in complexity, 
below we list one of the 3D EFIT normal stress equations, T n , for the case of anisotropic 
inhomogeneous material. Here we also point out that the ease of dealing with boundary 
conditions in EFIT makes it ideal for incorporating composite delaminations and other flaw 
types. 
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5 Conclusion 

We have developed a custom parallel 3D isotropic EFIT code to investigate Lamb wave 
scattering from damage in aerospace materials. We have found that 3D EFIT simulations can 
provide unique insight into complex Lamb wave behavior. The isotropic code implementation 
has been verified through comparisons with analytic solutions and has been compared against 
experimental results. In future work we plan to further develop and validate a 3D anisotropic 
EFIT code intended for modeling waves in composite materials. These computational tools 
will assist in the optimization of structural health monitoring flaw detection techniques for 
aerospace materials. 
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